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ABSTRACT 


A  method  is  developed,  using  the  quadratic  form  of  canonical  un¬ 
coupled  state  variables  as  a  Lyapunov  function  to  determine  the  switch- 
ing  logic  for  a  quasi-optimum  control  of  a  non-linear  dynamic  system 
Analytio  arguments  are  presented  to  support  the  method,  which  is  capable 
of  being  extended  to  any  order  system  and  any  number  of  controls,  sub- 
jeot  to  oertain  practical  limitations  noted  in  the  analysis . 

A  computational  scheme  for  determining  the  canonical  state  variables 
and  control  functions  is  presented  and  a  topological  interpretation  is 
made  of  the  choice  of  variables  used  for  the  control  logic  calculation.. 

The  controller  based  on  the  method  described  is  applied  to  a  non¬ 
linear  second-order  system.  Phase  trajectories  for  both  the  uncontrolled 
and  controlled  systems  are  obtained  by  means  of  a  digital  computer  simula¬ 
tion.  Various  aspects  of  the  theoretical  limitations  of  the  method  pre¬ 
sented  are  investigated  and  the  experimental  results  are  analyzed  with 
respect  to  the  theoretical  predictions. 

#The  system  is  postulated  to  be  described  by  a  system  of  differential 
equations  of  known  form. 
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lo  Introduction 


The  stability  of  dynamic  systems  whose  equations  of  motion  are  such 
that  their  solutions  are  difficult  or  impossible  to  obtain  in  closed  form 
may  successfully  be  analyzed  by  means  of  Lyapunov's  Second  Methods 

We  have  a  given  dynamic  system  described  by  the  set  of  differential 
equations  as 


x  *  i,(x);  jl(°)  s  0  °  (1-1) 

Then*  referring  to  Fig,  (l-l),  the  origins,  according  to  LaSalle  and 
Lefschetz®  is? 

stable  whenever  for  each  R  there  is  an  r  «fR  such 
that  if  a  path  (a  motion)  g*  initiates  at  a  point  x°  of 
the  spherical  region  S(r)  then  it  remains  in  the  spherical 
region  S(R)  ever  afters  that  is,  a  path  starting  in  S(r) 
never  reaches  the  bouncfeiy  sphere  H(R)  of  S(R)s 

asymptotically  stable  whenever  it  is  stable  and  ia  addition 
every  path  g^  starting  inside  some  S(RC),  R0^?  o,  tends  to  the 
origin  as  time  increases  indefinitely! 

unstable  whenever  for  some  R  and  any  r,  no  matter  how  small, 
there  is  always  in  the  spherical  region  S(r)  a  point  x 
such  that  the  path  g-*  through  x  reaches  the  boundary  sphere 
H(R)  e 


Figo  (l-l)5 


Geometric  Interpretation  of  the  Various  Types 
of  Stability, 
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Now  suppose  we  have  a  given  dynamic  system  described  by  the  set  of 


differential  equations  written  in  vector  form  as 

x  =  £(x,u) 


(1-2) 


or 

x  =  Fx  +  Du  +  £  (x,u)  (1-5) 

where  F  is  a  constant  square  matrix,  D  is  the  control  distribution 

* 

matrix,  ^  is  the  control  function  vector,  and  £  (x,u)  contains  terms  of 
a  second  and  higher  power  in  x  and  u.  The  linearized  form  of  the  equa- 
tion  is 

x  =  Ex  +  Du  (1-4) 

and  can  be  said  to  adequately  describe  the  motion  of  the  system  for 
small  values  of  x. 

According  to  Kalman  and  Bertram^,  a  Lyapunov  function  may  be  de- 

fined  as  some  scalar  function  of  the  state  variables  V(x)  with  the  fol¬ 
lowing  properties:  (a)  V(x)  )  0,  V(x)(  0  when  x  £  x0,  where  Xg  is  the 

o 

equilibrium  state,  and  (b)  V(x)  =  V(x)  s  0  when  x  -  xqC 

We  wish  to  apply  Lyapunov' s  Second  Method  to  the  design  of  a  con¬ 
troller  so  as  to  return  the  state  vector  (x)  to  a  position  of  equilib¬ 
rium  (Xg)  in  some  ootimum  manner.  No  loss  of  generality  will  result 
if  we  choose  (xg)  to  be  the  origin  of  the  state  space,  and  henceforth 
the  equilibrium  position  will  be  considered  as  such. 

Lyapunov  stated  that  if  a  Lyapunov  funotion  V(x)  exists  (being 
positive  definite)  in  some  neighborhood  of  the  origin,  and  if  -V(x) 
is  also  positive  definite  in  some  neighborhood  of  the  origin,  the 
system  is  asymptotically  stable^. 

In  general  the  method  employed  will  be  to  make  V(x)  as  negative  as 

possible  in  order  to  drive  the  system  to  the  origin  as  quickly  as 
possible. 

-  2 


The  Second  Method  of  Lyapunov  leaves  open  to  the  individual  inves¬ 
tigator  the  choice  of  a  particular  Lyapunov  function  to  be  used  in  the 
analysis  of  a  given  system*  For  a  given  system  then*  there  are  possibly 
an  unlimited  number  of  such  functionso  The  quadratic  form  of  the  canoni¬ 
cal  state  variables  as  the  particular  Lyapunov  function  employed  in 
the  design  of  the  controller  was  governed  by  the  simplification  obtained 
in  the  control  logic* 

Analysis  and  design  for  a  controller  of  a  linear ,  stationary  system 
using  the  quadratic  form  procedes  in  a  straightf orward  manner 0  In  the 
case  of  a  non-linear  system  however,  the  limitation  that  the  linear  ap¬ 
proximation  is  sufficiently  accurate  only  within  some  region  surrounding 
the  equilibrium  position  poses  a  rather  severe  restriction  on  the  extent 
of  the  state  space  for  which  a  given  control  can  be  useful* 

In  order  to  overcome  this  limitation  to  some  degree  the  controller 
was  designed  in  such  a  manner  that  the  control  logic  was  obtained  by 
successive  linear  approximations  to  the  actual  system  at  points  on  the 
system  trajectory# 

For  best  performance  of  the  system  the  time  involved  in  the  calcula¬ 
tion  of  the  control  logic  should  be  infinitesimal  so  that  there  would  be 
no  appreciable  time  delay  in  the  application  of  the  proper  control  func¬ 
tion*  That  is,  the  system  state  vector  would  not  have  moved  from  the 
point (x^),  where  the  state  variables  were  sampled,  to  a  new  position  on 
the  trajectory,  (x£)  where  the  control  function  was  actually  applied 0 

In  actual  practice  a  finite  calculation  time  is  of  course  necessary 
so  that  the  state  vector  will  always  have  moved  from  the  point  of  sampling 
before  the  new  value  of  control  function  can  be  applied# 

For  plants  with  long  time  constants,  such  as  might  exist  in  chemical 
or  metallurgical  processes  or  ship  stabilization  systems,  presently 
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available  digital  computers  have  computation  speeds  adequate  to  success*, 
fully  provide  on  line,  real  time  control. 

Application  of  this  method  to  plants  with  short  time  constants  will 
depend  upon  increased  speeds  of  the  computers  and  the  development  of  more 
refined  methods  of  calculation. 
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2.  Theory. 

If  we  have  given  a  Lyapunov  function  V(x)  s  x^Px  we  may  assert 
that^  the  equilibrium  state  xe  »  0  of  a  linear  dynamic  system 

x  *  Fbc  (2-1) 

is  asymptotically  stable  if  and  only  if  given  any  symmetric  positive- 
definite  matrix  Q,  there  exists  a  symmetric,  positive  definite  matrix  P 
which  is  the  unique  solution  of  the  set  of  linear  equations  described  by 

pTp  +  pp  =  — Q  .  (2-2) 


A  simple  procedure  for  calculating  the  control  logic  attributed  to 
R.W.  Bass^  procedes  as  follows:  Choose  Q^O  arbitrarily.  P>0  may  then 
be  calculated  from  (2-2).  ||x||2p  will  then  be  a  Lyapunov  function  for  the 

system  (2-1).  We  now  choose  u(t)  so  as  to  make  V  (x)  as  negative  as 
possible.  Let 

V(x)  »  x^Px  =  ||x||2P,  then 

V(x)  a  X^Px  +  x^Px  •  (2-3) 

We  wish  to  consider  the  effect  of  control  on  V.  Taking  the  transpose 
of  (1-4)  we  obtain 

xT  »  xTFT  +  u^D1  (2-4) 


and 


V(x)  ~  x^F^Px  +  x^PFx  +  u^D^Px  +  x^PDu 


or 

V(x)  =  xT(FTP  +  PF)x  4  2uTDTPx  . 

T 

Since  P  is  symmetric,  DP  *  PD.  Simplifying  this  gives 

V(x)  »  -x^Qx  +  2u_^D^Px 
or  V(x)  *  “||X||^Q  ♦  2u^D^ftc  . 

To  make  V(x)  as  negative  as  possible  we  want  the  term  2u  D  Px  to  be 
as  negative  as  possible.  Therefore  the  control  logic  is  given  as 


(2-5) 


5 


u^(t)  s  -a^  sgn  (D^Px)^  (2-6) 

where  the  are  constants  denoting  the  maximum  allowable  magnitudes  of 
each  individual  control  effort  u^.  It  is  obvious  from  (2-6)  that  in 
order  to  make  the  term  2ir-DAPx  as  negative  as  possible  the  individual 
controller  u^(t)  must  always  exert  a  maximum  effort ,  and  in  such  a 
direction  that  the  tern 

-aj_  sgn  (^Px)* 

is,  in  fact,  negative.  Thus  the  control  logic  becomes  merely  logic  for 
detecting  the  appropriate  switching  points  of  the  "bang-bang"  controller 
Ui(t). 

In  the  case  where  there  is  only  one  controller  the  vector  _u  may  be 
of  the  form  such  that  uj  *  (0,0,0  .  •  .  u).  The  system  equation  is 

then  of  the  form 


(2-7) 


~  0  1  0  ...  0 

0 

0  0  1 

0 

o 

o 

o 

*  + 

• 

• 

• 

1 

• 

-bn  -bn_i  .  .  .  -b^ 

U 

—  — 

_  _ 

where  the  b^  are  the  coefficients  of  the  characteristic  equation  of  the 
unforced  system,  taken  in  reverse  order. 

If  we  make  the  equivalence  transformation  from  the  physical  state 
space  to  a  canonical  state  space  denoted  by  the  matrix  equation 

Z  =  Ox  (2-8) 

then  by  differentiating  we  get 

Z  3  G*  (2-9) 

and  equation  (l-4)  becomes 

£  =  G^FGz  +  G_1Du  . 
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We  choose  G  such  that 


G_1FG  =  J\. 


(2-10) 


where  -A.  is  a  diagonal  matrix  with  the  same  eigenvalues  as 

For  the  system  with  one  controller  such  as  described  by  (2-7)  we 
may  further  choose  the  particular  transformation  matrix  G  such  that 

G*1Du  =  E 


where  s  (a,  a,  a,  .  .  .a).  This  has  the  effect  of  weighting  the 

value  of  all  the  state  variables  equally  in  determining  the  sign  of  u(t). 

If  the  P  of  (2-2)  is  chosen  to  be  the  identity  matrix  I  and  the  «A= 
matrix  substituted  for  F,  equation  (2-2)  becomes 

X T!  +  1  A.  «  -Q.  (2-11) 


If  all  the  real  parts  of  the  eigenvalues  of  F  are  negative,  so  then  are 
all  those  of  _A_ ,  and  Q  is  positive  definite.  Also 

$(x)  =  £T(_A_T  +  -A,  )%_  -fr  2ET£  .  (2-12) 

The  first  term  on  the  right  side  of  (2-12)  is  again  negative  and  the 
control  logic  is  given  by  the  second  term.  For  a  single  controller, 
we  have 

n 

u(t )  -  -a  sgn  (  £  y±)  •  (2-13) 

In  other  words  the  switching  function  is  such  that  the  control  must 
change  sign  whenever  the  sum  of  the  n  canonical  state  variables  changes 
sign. 


Consider  the  given  form 

“  -A.Z  +  g”1d^  • 


#It  is  always  possible  to  perform  the  equivalence  transformation  (2-10) 
provided  F  has  no  repeated  eigenvalues.  A  system  described  by  an  nth  order 
differential  equation  having  repeated  roots  may  be  transformed  into  a  par¬ 
tially  uncoupled  set  of  1st  order  differential  equations.  It  will  be 
assumed  that  the  systems  considered  in  this  paper  have  no  multiple  roots. 
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If  we  take  the  LaPlace  transform  we  get 

sY(s)  -  Z(0)  =  _A_Y(s)  +  G“^-D^(s) 

or 

Y(s)(sl  -  J^)  =  ^(0)  *  GT ^D_u( s )  .  (2-14) 

Since  from  (2-6)  the  control  variable  should  assume  only  values  of  +  a 

or  -  ^  ,  we  stipulate  that  it  is  a  constant  during  a  given  control 

period  (between  switchings ).  We  then  have  the  LaPlace  identity 

(  G“  lD_u)  ( s )  —  _Z_  • 

s 

Solving  (2-14)  for  a  given  component  y^  we  get 

Yi(s)(s  -Xj.)  -  yi(o)  +  — — 


or 


Yi(s)  = 


Yi(0) 


(s  -\±) 

Then  the  time  domain  solution  of  (2-15)  is 

y*(t)  -  yi(0)e*^it  - 


s(s  -X.) 


(2-15) 


*i 


(1  -eXit) 


or 

yj.(t)  =  -fyi(°)  +  --i-l  eX  1  -  —  .  (2-16) 

*i  J  X  ^ 

An  interesting  result  becomes  apparent  from  this  solution  as  t  in¬ 
creases  without  bound.  For  the  case  where  X^  is  negative  (a  consequence 

of  Q  being  positive  definite)  we  see  that 

,  .  a. 

yi(t)  -  -  }  t  oo 

xi 

This  result  indicates  that  although  this  form  of  controller  is  useful  in 
making  'O'(x)  more  negative  in  an  asymptotically  stable  system,  it  will  not 
result  in  V(x)  -*>  0  as  t  -►  OO  if  u  is  a  function  of  the  type  seen  in 


Fig.  (2-1). 
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/** 

U 
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Control  Function 

- 7-  l  U  J  1, 

As  An  Ideal  Relay. 

In  fact  it  will  result  in  a  condition  of  chatter-motion  around  the  origin 
as  the  controller  becomes  dominant  at  points  in  the  state  space  where  the 
yi<||e||s  where  ||e||  assumes  some  particular  small  value. 

This  objectionable  feature  may  be  overcome  in  the  case  of  asymptoti¬ 
cally  stable  systems  by  choosing  the  control  function  of  the  form  shown 
in  Fig.  (2-2)  below  so  that  the  controller  is  effectively  disengaged 
within  some  region  arbitrarily  close  to  the  origin,  allowing  the  system 
to  reach  the  equilibrium  condition  in  an  unforced  manner. 


Fig.  (2-2)  Control  Function  As  A  Relay  With  Dead  Zone. 

If  the  uncontrolled  system  were  such  that  one  (or  more)  roots  con¬ 
tained  positive  real  parts,  (in  this  case  -Q  can  no  longer  be  positive 
definite)  then  for  some  particular  points  in  state  space,  it  is  still 
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possible  that 


v(x)  =  -||x2||q  . 

However  V(x)>0  and  remains  so* 

Consider  a  plant  governed  by  the  equations 
dx/dt  =  Fx  +  Du(t) 

where  the  control  variables  are  subject  to  the  constraints, 

|ui(t)|  <  <  co  (i  s  1,  .  .  m)  .  (2=17) 

This  is  essentially  the  relay  or  saturating  servo  problem* 

The  problem  is  to  return  every  initial  state  to  the  origin* 

It  is  well  known  that  if  F  has  eigenvalues  with  posi¬ 
tive  real  Darts  then  there  are  some  states  which  cannot  be 
returned  to  the  origin  by  any  control  subject  to  the  con¬ 
straints  (2-17)4. 

If  a  controller  of  the  form  described  is  applied  to  the  unstable 
system,  it  is  apparent  from  (2-5)  that  control  may  be  maintained  so  long 
as  the  inequality  (2-18)  below  exists.2  That  is 

-||x||2Q  +  2utDtPx  <  0.  (2-18) 

From  the  solutions  (2-16)  the  behavior  of  each  y^(t)  for  0<t<k, 
where  k  denotes  the  switching  period  is  as  in  Fig*  (2-3). 


Fig.  (2-3)  Plot  of  y^(t)  vs.  Time  For  u^  =  -a^  sgn  (jf^i* 

Stable  Eigenvalue. 
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Fig.  (2-4)  Plot  of  y^(t)  vs.  Time  For  *  -a^  sgn  (E^)^ 

Unstable  Eigenvalue,  y-;  (0)  >  |  ai  |. 

I  Xi  I 


m 

Fig.  (2-5)  Plot  of  y^(t)  vs.  Time  For  =  -ai  sgn  (E.^)^ 
Unstable  Eigenvalue,  yj(0)  <|  ai  1  . 

I  X  i  I 

In  the  event  that  the  system  distribution  matrix  D  were  such  that 
the  controls  were  uncoupled,  i.e.  for  each  y^  there  were  a  corresponding 
u^,  then  it  would  be  possible  to  design  a  controller  whose  individual 
elements  changed  sign  at  the  zero  crossings  of  the  individual  state 
variables. 
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From  Figo  (2-5)  it  can  be  seen  that  the  element  would  switch  at 
t  s  To  It  is  possible  to  calculate  *  ^(y-j^O))  from  (2-16)  for 

the  case  where  y^(0)>0#  <0  . 

yi(t)  s  J  yi(0)  -  _ !^L=le^iT  *  — ii-  ,  let  yA(t)  s  0 

—  i  Jr  i 

then 


giving 


Xi 

■  ■  iT  —  ■  ■  —  — . 


and 


(2-19) 


From  (2-19)  it  is  clear  that  for  y^(0)  *  0*  T  s  — _ — .  ln(l)  -  0 

indicating  that  this  system  also  will  result  in  chatter-motion  around 
the  origin. 

Selection  of  a  control  function  of  the  form  of  Fig,  (2-2)  would 
be  unsatisfactory  in  this  case  because  the  uncontrolled  system  is  un¬ 
stable  at  the  origin,  A  better  form  of  function  would  be  that  of  the 
form  of  a  saturating  amplifier  with  a  high  gain0  (See  Fig,  (2-6)). 

Aui  ' 

-  (ins). 

A> 

”  -  4^ 

! 

Fig,  (2-6)  Control  Function  With  the  Characteristics  of  a 
Saturating  Amplifier. 
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Employing  LaSalle  and  Lefschetz's  definitions^  of  stability  des- 
cribed  in  Section  1»,  a  system  whose  eigenvalues  all  have  negative  real 
parts  is  asymptotically  stable  with  no  controls  merely  stable  with  a 
"bang-bang"  controller  (Figo  (2-l)),  and  again  asymptotically  stable  if 
the  controller  contains  a  dead  zone  as  in  Fig#  (2-2) o  Given  an  arbit¬ 
rary  H(R),  a  system  with  eigenvalues  lying  in  the  right  half  plane  is 
defined  as  unstable,  and  the  same  system  with  a  controller  of  the  form 
described  may  be  defined  as  stable  if  X°  is  restricted  to  lie  within 
some  S(r)  whose  radius  r  is  of  some  necessarily  small  value#  The 
restrictions  on  the  allowable  values  of  r  are  necessarily  related  to 
•the  constraints  on  the  magnitudes  of  the  control  functions  (2-17)» 

Given  a  non-linear  system  described  by  the  equation  x  *  F(x)x, 
the  linear  approximation  of  the  non-linear  system  referred  to  by  Kalman 
and  Bertram  may  be  obtained  by  a  Taylor  series  expansion  of  the  coef¬ 
ficient  matrix  F(x)  at  the  origin#  The  linearized,  constant,  matrix  F 

^  P „  / 

has  the  elements  — i - =  «  (The  matrix  F  is  the  Jacobian  of  F(x) 

o  Xf 

with  respect  to  x  at  x  s  x(0)  *  0)  « 

Suppose  a  system  described  by  the  equation 
x  =  F(x)x  +  Du 

.is  to  be  controlled  in  an  environment  where  the  disturbances  are  of 
such  a  magnitude  as  to  carry  the  state  variable  beyond  the  region  of 
validity  of  the  linear  approximation  given  by  (l-4)#  Referring  to 
Fig#  (2-7),  this  condition  corresponds  to  an  x°  lying  outside  of  S(r)0 
By  expanding  F(x)  about  the  point  x°  and  obtaining  the  Jacobian  of 
F(x)  with  respect  to  x  at  x  s  x°  we  may  obtain  a  linear  approxima¬ 
tion  at  x° 


c  ^ 

x  s  F  x  ♦  Du 


(2-20) 
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Figo  (2-7)  Geometric  Interpretation,  of  Successive  Linear 
Approximation  Process « 

Control  logic  based  on  this  approximation  will  be  valid  while  the 
trajectory  remains  within  the  spherical  region  S(p)<»  If  S(pn)  is  the 
region  of  validity  of  the  linear  approximation  due  to  the  n^*1  sampling p 
then  by  an  iterative  process  such  that  the  trajectory  from  Xa  to  Xn+“ 
always  remains  inside  S(pa)  we  may  obtain  an  optimum  control  policy 
based  on  the  particular  Lyapunov  function  employed  for  the  complete 
trajectory  for  x°  to  xe  *  0o 
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3 e  Application  of  Theory  to  Design  of  a  Single  Controller* 

It  was  shown  in  Section  2*  that  the  control  logic  resulting  from 

the  choice  of  the  quadratic  form  of  the  uncoupled  state  variables  as 

the  system  Lyapunov  function  is  of  the  form 

u(t )  «  -a  sgn  |  i^p1  y* 

when  the  system  has  only  one  controller®  Thus,  the  control  problem  is 

reduced  to  determining  the  sign  of  the  sum  of  canonical  state  variables 

and  insuring  that  the  control,  u(t),  assumes  the  opposite  sign..  The 

steps  involved  in  this  process  consist  of  the  followings 

I®  Measuring  the  physical  state  variables  (x)« 

2«  Calculating  the  canonical  state  variables  {%)  from  the 
measured  (x). 

3.  Forming  ^y^  • 

4.  Applying  the  proper  control  u(t)  ^  -a  sgn^yj.® 

As  was  shown  in  Section  2.,  the  variables  (^)  may  be  calculated  by 
the  equation 

£  -•  G^x  .  (3-1) 

-1 

The  problem  then  reduces  to  evaluating  G  , 

It  may  be  shown  that  a  matrix  G  such  that 

G"1FG  =  _A_  (3-2) 

v<here_/\-is  a  diagonal  matrix,  may  be  formed  by  calculating  the  eigen¬ 
vectors  of  the  matrix  F  and  constructing  an  array  in  which  each  eigen¬ 
vector  is  a  particular  column  of  the  array#  Once  a  G  has  been  obtained 
it  is  merely  necessary  to  calculate  the  inverse,  G~~,  in  order  to  be  able 
to  calculate  the  state  variables  (jr)# 

The  elements  of  the  diagonal  matrix  j\=  are  the  eigenvalues  of  F# 

If  the  system  is  oscillatory  the  elements  of  _A. will  necessarily  be 
complex.  Since  the  elements  of  the  coefficient  matrix  F  of  the  original 
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differential  equation  are  all  real,  consideration  of  (3=2)  reveals  that 
for  an  oscillatory  system  there  must  be  complex  elements  in  G  and  G“^s 

When  G  is  complex,  care  must  be  taken  in  evaluating  G“^t  If  we 
let  G  a  (feG  +  j$G,  and  attempt  to  form  G*1  -  (&G)-1  ♦  j(<3.  G)*1 

difficulty  arises  immediately0  Since  complex  roots  arise  in  conjugate 
pairs,  there  will  be  two  eigenvectors  (columns)  of  (t&G)  due  to  the 
identical  real  parts  of  the  conjugate  pairso  These  vectors  will  be 
identical,  and  therefore  linearly  related  so  that  (CfeG)  will  be  singulars 
Similar  reasoning  shows  that  (c$  G)  will  also  be  singular,  so  that  it  is 
impossible  to  obtain  (&G)“'^  and  (<D  G)~^a 

Let  us  form  a  new  matrix  H  in  the  f ollovdng  manner  e  If  the  order 
of  G  is  n,  the  order  of  H  will  be  2n  so  that  for  every  complex  element 
of  G  there  will  be  four  elements  of  H  arranged  in  a  square  pattern©  The 
two  elements  on  the  principal  diagonal  of  the  square  have  a  value  equal 
to  the  real  part  of  the  corresponding  complex  element  of  Go  The  ab¬ 
solute  values  of  the  remaining  two  elements  are  equal  to  the  imaginary 
part  of  the  corresponding  complex  element  of  G,  the  element  in  the  lower 
left  hand  corner,  taking  the  positive  sign,  while  that  in  the  upper  right 
hand  corner,  the  negative  signe  Where  the  corresponding  element  in  G 
is  real,  zeros  appear  in  the  locations  in  H  corresponding  to  the  missing 
imaginary  parts 0 

If  now  we  form  then  we  may  obtain  the  elements  of  from 

the  proper  locations  of 

However  since  we  are  interested  in  obtaining  the  (y)  corresponding 
to  the  measured  (x),  the  simplest  procedure  is  to  form  the  2  x  a  sup- 
plementary  vector  (x)  which  has  alternately  the  real  and  imaginary 
parts  of  the  physical  state  variableso  Since  these  state  variables  may 


IS 


be  measured  and  have  physical  significance,  they  may  only  be  real,  so 

♦ 

that  (x)  consists  of  the  values  of  the  state  variables  alternating  with 
zeros.  Matrix  multiplication  of  H”4  on  the  right  by  ( x )  results  in  a 
2  x  n  column  vector  whose  elements  are  the  real  and  imaginary  parts  of 
the  canonical  space  variables  (%) , 

EXAMPLE,  The  product  of  a  complex  2x2  square  matrix  and  a 
column  vector  with  real  elements 

R,,  +  .iln  Ri 


n 

yz 


Then 


*1 
yz 

Rearranging  gives 


‘11  4  j^ll  R12  4  3j12 

R21  +  JI21  **22  4  ^22 


Rllxl  4  ^ll3^  4  R12x2  4  JI12X2 

R21x1  +  jI21Xl  4  R22X2  4  jI22X2 


x2 


yl  - 

(Rnxi  + 

r12x2) 

4  j(ln^i  +  Ii2x2^ 

(3-3) 

y2  3 

4 

^ZZXZ^ 

4  *j^I21Xl  4  J22x2) 

(3-4) 

Suppose  instead  we  form  a  4  x  4  matrix  of  the  form  H"*,  Then 


&yi 

Rn 

*'jIll 

R12 

-^12 

X1 

&Yi 

4^n 

R11 

r12 

X 

0 

&y2 

R21 

^’I21 

R22 

-jl22 

x2 

<*y2 

3^21 

% 

*^I22 

R22 

L°. 

and 

ftyi  3  Rnxi  4  hzxz 

Ayi  s  Pn*i.'  4 

ay2  s  ®21xl  4  R22x2 

*  A  y2  ■  JI21X1  +  jl22x2  • 

Rearranging  gives 

yi  s  (Rllxl  4  R12x2^  4  j^llxl  4  I12x2^ 
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V2  *  (R21xl  +  R22x2^  +  ^I21xl  +  I223C2^ 

which  is  identical  to  (3-3)  and  (3-4). 

It  is  evident  from  the  above  discussion  that  in  the  case  of  an  oscil¬ 
latory  system  the  transformation 

%  *  G_1£  (3-5) 

may  map  a  given  physioal  state  variable  x^  from  a  point  on  the  real 
axis  to  some  point  in  the  complex  plane.  Thus  for  a  given  coordinate 
axis  in  the  physical  n-space  there  is  a  complex  plane  in  the  canonical 
spaoe. 

Since,  however,  there  is  a  one  for  one  correspondence  between  points 
in  the  two  spaces,  a  control  which  succeeds  in  reducing  all  state  vari¬ 
ables  to  zero  in  the  canonical  space  will  have  also  reduced  the  real 
state  variables  to  zero. 

The  computation  method  employed  to  determine  the  transformation 
matrix  G  was  based  on  a  computer  program  MA.TSUB,  programmed  by 
Louis  W.  Ehrlich  of  the  University  of  Texas,  which  calculates  eigen¬ 
values  and  eigenvectors  of  a  matrix  up  to  order  50.  The  general  method 
employed  is  to  make  an  initial  guess  for  each  eigenvector  of  the  matrix 
P  and  to  converge  on  the  correct  value  by  an  iterative  scheme  due  to 
E.E.  Osborne  .  The  original  program  employed  a  guess  of  1.0  for  all  com¬ 
ponents  of  the  eigenvectors.  Since  it  was  assumed  that  the  system 
matrix  F  was  not  changing  at  more  than  a  moderate  rate,  a  modification 
was  made  to  the  program  so  that  the  most  recently  calculated  value  for 
each  eigenvector  was  selected  as  the  initial  guess  for  a  new  calculation 
thereby  speeding  the  computation.  The  output  of  the  MATSUB  program  was 
also  modified  so  that  the  eigenvectors  were  arranged  into  two  n  x  n 
matrices  corresponding  to  the  real  and  imaginary  parts  of  the  transfor¬ 


mation  matrix  G 
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From  the  two  n  x  n  matrices  the  2n  x  2n  matrix  H  was  formed  and  the 


inverse  taken.  The  canonical  state  variables  were  then  determined  and 
the  proper  control  calculated  and  applied. 

The  path  of  the  system  trajectory  in  the  physical  space  (x)  was 
then  calculated  for  a  selected  time  interval  between  sampling  instants 
by  the  Runge-Kutta-Gill  method.  Instantaneous  values  of  the  elements 
of  the  coefficient  matrix  F  =  F(x)  were  calculated  for  the  new  sam¬ 
pling  point  and  the  procedure  repeated  until  the  trajectory  in  the 
canonical  space  was  within  a  pre-set  epsilon  region  of  the  origin  or 
a  given  time  had  elapsed. 

Initially  the  ideal  case  was  simulated  by  holding  the  solution  of 
the  system  trajectory  at  the  sampling  point  until  the  new  calculation 
of  the  control  function  was  made  and  the  control  applied.  Simulation  of 
the  non-ideal  case,  a  finite  calculation  time,  was  affected  by  applying 
the  control  calculated  at  the  (n-l)^  sampling,  at  the  time  of  the  n^b 
sampling.  This  is  equivalent  to  a  controller  continually  applying  a 
new  control  as  soon  as  it  can  perform  the  sampling  and  control  cal¬ 
culations  . 

Simplified  flow  charts  of  both  the  ideal  and  non-ideal  case  simula¬ 
tions  are  presented  in  Appendix  I.  The  computer  program  in  FORTRAN 
language  (for  the  ideal  case)  is  also  included  in  Appendix  I. 
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4.  Simulation  of  Control  of  a  Non-Linear  Second  Order  System. 

A.  Background. 

The  non-linear  VanderPol  equation 

X  -  (1  -x2)  x  +  x  *  0  (4-1) 

was  selected  as  the  system  on  which  to  evaluate  the  controller.  Al¬ 
though  it  is  only  a  second  order  equation  it  was  felt  that  since  it 
is  quite  non-linear  it  would  be  a  fair  test  of  the  controller  and  the 
two  dimension  phase  space  has  the  advantage  of  making  graphical  analysis 
more  simple. 

The  equation  exhibits  a  stable  limit  cycle  of  the  general  shape 
shown  in  Fig.  (4-1).  • 


Fig.  (4-1)  Phase  Portrait  For  an  Uncontrolled  VanderPol  Equation. 
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; 


Referring  to  (4-l),  when  x  ■  0,  the  system  is  most  unstable.  For 
this  point  the  roots  are  +.5  ±  j.866.  It  was  determined  from  trajec¬ 
tories  of  the  uncontrolled  system  obtained  by  digital  computer  analysis 
that  when  the  system  was  in  the  limit  cycle  the  maximum  excursion  of 
displacement  was  ~  2.1925.  This  value  was  substantiated  by  an  analog 

computer  simulation.  From  (4-l)  the  eigenvalues,  or  roots,  of  the  system 
vrhen  the  displacement  reaches  its  maximum  value  in  the  limit  cycle  were 
calculated  to  be  -3.5235  and  -.2835  so  that  the  locus  of  roots  of  equa¬ 
tion  (4-1 )  during  a  complete  circuit  of  its  stable  limit  cycle  is  shown 
in  Fig.  (4-3).  .  • 


Fig.  (4-3)  Root  Locus  of  the  Uncontrolled  VanderPol  Equation 
in  Stable  Limit  Cycle. 

Examination  of  Fig.  (4-l)  or  equation  (4-l)  reveals  that  the  system  be¬ 
comes  unstable  whenever  |x|  <  1.0  and  becomes  stable  again  whenever 
|x|>  1.0.  Thus  the  equation  presents  the  controller  with  the  problem 
of  controlling  a  system  which  contains  at  various  times,  stable  real 
roots,  stable  complex  roots,  and  unstable  complex  roots. 

It  is  informative  to  observe  that  a  stability  analysis  of  the 
VanderPol  equation  by  Lyapunov's  Second  Method  yields  precisely  the 
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same  results  as  the  more  familiar  methods  of  analysis  such  as  the  root 
locus. 


Rewrite  (4-1)  in  vector  form.  Then 


or 


and 


Choose 


then 


i  •  [-?  <^>]  z. 


*1  1  *2 


X2  =  +  (l-x^)x. 


2  * 


V(x)  =  j|x||  -  xj  +  *§ 


V(x)  ■  2x^xi  +  2x2xg 
and  from  (4-3)  and  (4-4), 

'O'(x)  =  2x^x2  +  2x2(-x^  +  (1-Xj)x2) 


(4-2) 

(4-3) 

(4-4) 


(4-5) 


or 

V(x)  =  2x| ( 1-x^ )  (4-6) 

so  that  V(x)  is  positive  definite  but  -V(x)  is  not  for  |x^j<l.  There¬ 
fore  the  region  near  the  origin  (-l<x<l)  is  unstable.  This  is  pre¬ 
cisely  the  result  obtained  by  observing  that  the  root  locus  passes  into 
the  right  half  plane  whenever  |xjJ  <1. 

B.  Procedure. 

It  was  decided  to  initially  investigate  the  action  of  the  control 
with  no  time  delay  between  calculation  of  the  control  and  its  application 
in  order  to  verify  that  the  method  of  control  and  the  simulation  pro¬ 
cedure  (described  in  Section  3.)  were  feasible.  First  a  family  of  tra¬ 
jectories  of  the  system  with  no  control  effort  was  obtained  from 
various  initial  conditions  chosen  inside  and  outside  the  limit  cycle 
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(Graphs  A-l  through  A-8).  Trajectories  were  then  obtained  using  the  same 
initial  conditions  with  the  control  system  activated  (Graphs  B-l  through 
B-8).  The  sampling  rate  was  set  at  .3  seconds  and  u  was  set  at  1.0*  An 
arbitrary  epsilon  region  surrounding  the  origin  was  established  so  that 
the  trajectory  was  considered  to  have  reached  the  origin  when  the  sum 
of  the  canonical  state  variables  squared  was  less  than  .001.  Whenever 
the  trajectory  entered  the  epsilon  region  of  the  origin  the  solution 
was  terminated. 

Effects  on  the  trajectories  due  to  varying  sampling  rates  (Graphs 
C-l  through  C-5)  were  then  investigated,  and  finally  two  trajectories  were 
obtained  for  the  system  simulating  the  effect  of  finite  calculation  time 
(Graphs  D-l  and  D-2).  These  trajectories  all  were  started  from  an 
initial  condition  close  to  the  uncontrolled  limit  cycle  trajectory. 

The  simulation  program  included  a  sub-rort ine  for  the  graphical  pre¬ 
sentation  of  the  system  trajectories.  The  graphs  noted  above  are  pre¬ 
sented  in  Appendix  II. 
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UNCONTROLLED  CASE:  SAMPLING  RATS  s  0.3  SECOND. 
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II-CU  CONTROLLED  CASE:  SAMPLING  RATE  a  .50  SECONDS,  u  =  1*0 
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TIME-DELAY  CASE:  INITIAL  CONDITIONS  x  '  =  1.5,  X  =  1.5,  u  a  1,0 
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D.  Discussion  of  Digital  Computer  Simulation, 

Comparison  of  the  trajectories  of  the  controlled  and  uncontrolled 
systems  starting  from  the  same  initial  conditions  (Graphs  A-l  through 
A-8  and  B-l  through  B-8)  clearly  indicates  that  the  control  system  with 
a  control  effort  =  1.0  and  a  sampling  rate  of  .30  seconds^'  is  capable 
of  driving  the  system  into  the  unstable  region  surrounding  the  origin  and 
maintaining  it  within  some  region  of  lesser  extent  than  that  enclosed  by 
the  system’s  uncontrolled  limit  cycle.  Referring  to  Pig.  (l-l),  if  the 
radius  of  S(R)  were  established  for  this  system  as  R  *  .5,  then  the  con¬ 
trolled  system  might  be  defined  as  stable  while  the  uncontrolled  system 
would  be  defined  as  unstable.  Graph  B-l,  initial  conditions  x  ■  .2, 
x  »  .2,  illustrates  the  system  trajectory  in  the  region  close  to  the 
origin.  The  chatter-motion  mode  of  operation  is  readily  apparent. 

The  theoretical  analysis  of  Section  2.,  predicting  that  by  employ¬ 
ment  of  this  type  of  controller,  an  unstable  system  could  be  made  stable, 
but  not  asymptotically  stable,  and  that  a  chatter -mot ion  mode  of 
operation  would  result  around  the  origin,  is  substantiated. 

Examination  of  Graphs  C-l  through  C-5  reveals  that  an  increase  in  the 
control  system  sampling  rate  causes  the  trajectory,  having  once  arrived  in 
the  neighborhood  of  the.  origin,  to  remain  thereafter  within  a  smaller  re¬ 
gion  of  the  origin  than  the  same  system  operating  with  a  low  sampling  rate. 

Curves  B-3  and  C-3,  both  trajectories  for  the  identical  initial  condi¬ 
tions  and  system  parameters,  demonstrate  remarkably  different  trajectories, 
each  of  which  eventually  arrives  in  a  neighborhood  close  to  the  origin. 

#The  time  interval  between  successive  calculations  of  the  control 
function.  The  Runge-Kutta  integration  of  the  trajectory  employed  an 
interval  of  0.03  seconds. 
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Examination  of  the  calculated  data  revealed  that  in  the  case  of 


Graph  B-3  the  initial  values  of  the  canonical  state  variables  were 

=  .75  -  j.1561 

y2  =  .75  +  j.1561 

while  in  the  case  of  Graph  C-3  the  values  were 

y-L  ■  .75  -j.1561 

Y2  =  .75  —j.1561  « 

Since  the  control  function  u  *  -1.0  sgn  (  E&yi  +  E<!yib 

it  is  evident  that  u  assumed  a  different  sign  in  each  case.  That  this 
was  so  can  be  determined  by  examination  of  the  diverging  paths  of  the 
trajectories  in  the  two  cases. 

It  seemed  evident  that  the  more  direct  trajectory  should  be  the 
correct  one,  however  the  possibility  was  raised  that  there  might  not 
be  a  unique  solution  for  the  control  function.  Six  additional  runs 
each  exactly  coincided  with  the  more  direct  trajectory  B-3,  indicating 
that  oerhaps  an  error  in  calculation  led  to  the  selection  of  the  op¬ 
posite  value  of  control  function  rather  than  there  being  two  possible 
solutions  for  the  control  function. 

The  computation  routine  had  been  programmed  so  that  the  inverse 
transformation  matrix  was  printed  at  each  sampling  point.  It  was  dis¬ 
covered  that  the  transformation  matrices  in  the  two  cases  were 
not  identical. 


INITIAL  INVERSE  TRANSFORMATION  MATRIX  FOR  GRAPH  B-3 


.35820E-10 
-.64051E-00 
-  .35820E-10 
.640516-00 


.64051E-00 
.35820E-10 
-.64051E-00 
- .35820E-10 


. 50000E-00 
-.40032E-00 
.50000E-00 
.40032E-00 


.40032E-00 
. 50000E-00 
-.40032E-00 
.50000E-00 
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INITIAL  INVERSE  TRANSFORMATION  MATRIX  FOR  GRAPH  C-3 


.37180E-08  .64051E-00  .50000E-00  .40032E-00 

-.64051E-00  .37180E-08  -.40032E-00  .50000E-00 

©  50000E-00  .40032E-00  .65670E-09  .64051S-00 

- .40032E-00  .50000E-00  -.64051E-00  .65484E-09 

In  order  to  determine  if  two  solutions  were  possible,  the  two  in¬ 
verse  matrices  were  checked  against  the  original  matrix  to  see  if  both 
were  valid o 

The  original  matrix  was  calculated  by  hand  in  the  same  manner  that 
the  computer  was  programmed  to  do0 

The  procedure  is  illustrated  below. 


The  initial  F  matrix  was 


fo  1  1 

L-l  -1.25J 


and  the  initial  eigenvalues  were 


-.625  *  j. 7806.  Then  for  the  first  eigenvalue 

s  (-.625  *  o„7806) 


p  [:;] 


or 


x2 

(-.625 

*  j.7806)x^ 

*1 

* 

(-.625 

-  o»7806)x2 

so  that 

x2 

1.0 

X1 

= 

-.625 

-  o»7806  o 

In  a  similar 

manner,  for 

the  second 

eigenvalue 

x2 

S3 

1.0 

X1 

S 

-.625 

^  0  • 7 606  , 

Then  the  2n  : 

k  2n  matrix  becomes 

-.625 

-.7806 

-.625  .7806 

.7806 

-.625 

-.7806  -.625 

1 

0 

1  0 

0 

1 

0  1 

When  this  matrix  is  multiplied  by  the  inverses,  B-3  and  C-3,  only  the 
inverse  B-3  results  in  the  identity  matrix,  demonstrating  that  the 
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calculation  resulting  in  C-3  is  not  an  independent  correct  solution* 

Comparison  of  the  false  trajectory  with  the  correct  one  indicated 
the  possibility  of  there  being  at  least  two  other  instances  of  erroneous 
calculations  besides  the  one  at  the  origin*  These  are  indicated  by  the 
abrupt  changes  in  direction  (outward)  of  the  trajectory  at  the  points 
x  -  .65,  y  s  0,6,  and  x  £  '.5,  y  s  *45  .  Calculations  similar 
to  the  one  carried  out  above  verified  the  fact  that  erroneous  determina¬ 
tions  of  the  correct  control  function  had  been  made* 

The  similarity  of  the  trajectories  of  curves  C-l,  C-2  and  D«X  in 
the  fourth  quadrant  raises  the  interesting  comparison  between  the  control 
obtained  by  Lyapunov's  method  and  the  method  of  optimum  control  employing 
the  calculation  of  the  optimum  switching  lines  in  negative  time  as  dis¬ 
cussed  by  I.  Flugge-Lotz  and  H.A.  Titus^. 

We  may  speculate  that  there  is  some  optimum  switching  line  for 
this  system,  emanating  from  the  origin  in  a  manner  similar  to  that 
in  Fig*  (4-4).  * 


Fig.  (4-4)  System  Optimum  Switching  Line. 
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The  optimum  trajectory  for  this  system  would  then  proceed  from  the 
initial  condition  to  its  first  intersection  with  the  optimum  switching 
line,  change  the  sign  of  its  controller,  and  proceed  along  the  optimum 
switching  line  to  the  origin. 

It  appears  evident  that  the  method  of  switching  logic  discussed  in 
this  paper  selects  a  switching  time  too  soon  in  the  cases  of  curves  C~1 
and  C-2,  In  the  case  of  curve  D-2  the  system  appeared  to  be  in  a  state 
of  indecision.  The  switchings  at  points  A  and  C  undoubtedly  would  have 
brought  the  trajectory  initially  closer  to  the  origin  than  the  one 
finally  selected  at  E.  The  switchings  at  B  and  D  served  only  to  hinder 
the  system  response.  Investigation  of  the  computer  data  again  indicated 
erroneous  calculations  made  for  switchings  at  points  B,  C  and  D,  although 
the  resulting  decision  at  C  proved  to  be  a  correct  one. 

Graphs  D-l  and  D-2,  showing  the  simulation  of  the  system  operating 
under  non-ideal  conditions,  again  substantiate  the  fact  that  higher 
sampling  rates  generally  produce  more  satisfactory  operation  of  the 
system.  It  is  interesting  to  note  that  the  controlled  system  operating 
with  a  time  delay  of  .3  seconds  eventually  settles  into  a  stable  limit 
cycle  of  approximately  one  half  the  size  of  the  uncontrolled  system. 

The  effect  of  the  controller  on  the  excursions  of  the  system  roots 
as  the  state  point  travels  along  a  trajectory  for  a  particular  initial 
condition  is  demonstrated  in  Figs.  (4-5)  and  (4-6).  The  numbers  indicate 
time  in  seconds  to  reach  the  position  indicated. 

The  root  locus  for  the  uncontrolled  system  starting  from  initial 
conditions  x  =  .5,  x  =  .5  is  illustrated  in  Fig.  (4-5).  The  roots 

attain  their  most  stable  positions  after  the  trajectory  has  settled  into 
its  stable  limit  cycle. 
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Fig,  (4-5)  Root  Locus  for  Uncontrolled  System,  Initial 
Conditions  x  =  ,5,  x  ~  ,5, 

As  the  controlled  trajectory  approaches  the  region  of  chatter- 
motion  surrounding  the  origin,  the  root  locations  depart  very  little 
from  their  most  unstable  position  (+.5  ±  j«866)  because  maximum 

instability  occurs  when  x  ~  0,  The  most  stable  root  positions  on 

the  controlled  system  root  locus  occur  early  in  the  solution  while 
the  trajectory  is  relatively  far  from  the  origin. 


Fig,  (4-6)  Root  Locus  for  Controlled  System,  Initial 
Conditions  x  =  .5,  x  s  ,5. 


35 


5.  Conclusions* 


The  operation  of  the  controller  in  regulating  a  simulated  non* 
linear  system  described  by  the  VanderPol  equation  supports  the  theo¬ 
retical  predictions  developed  in  the  design*  The  solution  obtained  by 
this  method  is  seen  to  be  less  than  optimum  because  the  logic  does  not 
cause  the  control  to  switch  at  the  optimum  svdtching  curve*  It  is 
apparent  from  the  figures  that  the  sampling  rate  for  determining  the 
eigenvalues,  and  hence  changes  in  control  action,  is  critical  for  this 
to  be  an  effective  control  design  procedure*  With  the  fine  sampling 
the  system  disturbances  were  rapidly  and  effectively  controlled* 

There  exist  today  very  few  direct  techniques  for  the  control  of 
non-linear  systems*  The  procedures  studied  here  provide  such  a  tech¬ 
nique*  The  method  has  the  advantage  of  being  applicable  to  non* 
linear  time  varying  systems  of  any  order*  It  does  however  require  a 
digital  computer  in  the  control  link* 
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APPENDIX  I 

FLOW  CHARTS  AND  FORTRAN  LANGUAGE  PROGRAM 
FOR  DIGITAL  COMPUTER  SIMULATION 


37 


Fig.  1-1  FLOW  CHART  FOR  IDEAL  CASE  SIMULATION. 
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Fig.  1-2  FLOW  CHART  FOR  FINITE  CALCULATION  TIME  SIMULATION 
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TOEL  =  7-«-  DELTA 

1304  IF (ID EL  -  T  -  DT)  1300  »  1301  *  1301 

m  isibsv.iev302'™-*7* 

GRlT0’l000  ^  *  {U°  "  XRM)  *  XR (  1 )  )  i 

1303  CONTINUE  ■' 

PRINT  21202 

21202  FORMAT^  26h-{RAJ£CTQRY  INSIDE  EPSILON///) 

<  f\  l  W  I  n  J  Z  U  i! 

h,2°2  T  ;  YSCR  R4DIUS  VKT0R///’ 

>  777  CONTINUE  ’ 

GRAPH'fNUNPTS.XI.YI.B! 

C1S01UpLl'Riur|ArfwiC®NTR?Lx|0,SDfTfMU0P*JECT0RY  «*  «««■»«•  P£RI0D 

nos  1**05,  1406, 1406 

LINE  =  4 

1406  PRINT  1 202 , T, XR ( 1 ) , XR { 2 ) 

LINE  =  LINE  +  1 
I  Ft  900  -  NUXPT  S  5  1304,750,750 
7b0  NUMPTS  =  NUN? T S  +  1  ’ 

XH  VJMPTS)  -  X R  {  1  ) 

YI f  NUMPTS)  =  XR { 2 ) 

CO  f 0  1304 
END 

SUBROUTINE  RKUT7A  (  M,  T  ,  XR  ,DT  ,  UOP) 

CCNSopl1  ,£?*ggi&f?iAt6ll'xt>OT«2c>»  xcl20,t  Cl4’ 

C  {  1  )  —  0  o  0 

c ;  2  >  =0.5 

C  t  3 )  =  0.5 
C  (  4  1  =1.0 
DO  1050  I  =1,4 
TC  =  T  +  C  (  I  )  *  OT 
DO  1051  J  =  1  ?  M 

i051  XC(J)  =  XR  (  J  )  •}•  C  (  I  )  *  AK  {  I  —  1  ,j) 


CALL  DERI  V  ( TC  ,  XC  ,  XDGT 


UOP 


) 


41 


+  2.0*  AK(2,J)  +  2 . C* AK! 3  »  J  5  +AK(4,J))/6. 


+  A  R  (  M ,  2  )  *  XR  (  2  )  +  UOP 
ALIS  ) 


DO  1050  J  =UM 
1050  AX  (I.J)  =  DT  *  XDOT ( J ) 

DO  1052  J  =  1  »H 
1052  XR ( J )  =XR ( J )  +  ( AK ( 1 »  J ) 

END 

SUBROUTINE  DERI  V  (T,XR,  XDOT,  UOP,  M  ) 

DIMENSION  XDOT  (20),  XRI2C),  AR ( 20, 2 0 ) , A I (20, 20 ) ,GR ( 20 , 20 ) , G I ( 20 , 
10), VALUES! 2, 20) 

COMMON  AR  ,A1 ,GR,GI , VALUES 
XDOT! 1 )  =  XR ( 2 ) 

XDOT! 2)  =  AR ! M, 1 )  *  XR(l) 

END 

SUBROUTINE  GSUB  (M  »  ALRS  ,  AHi  i 
ODIMLNSIQN  AR(20,20) ,AI (20,20) tBR(20,  20) , B I ( 2 0 * 2 0 ) , CR! 20 , 20 ) , 

1CI !  20,20 ) »XR( 20 ) ,XI (20  5  , YR ( 20) , YI ( 20 ) , VALUES ( 2 , 20 ) , GR ( 20 , 20 ) , 

2  G I !  2  0 , 2  0  ) 

COMMON  CR  , Cl  ,GR»GI .VALUES 
N  =  M 
LOOP  =  1 

922  DO  519  I  =  1 » N 
DO  519  J  =  1 , N 
BR!  I  ,  J)  =  CR! I, J  ) 

AR!  I  - J )  =  CR!  I , J ) 

=  Cl  !  I  ,  J) 

519  A I ( I , J )  =  Cl  !  I, J  ) 

ASSIGN  917  TO  IA 
ID 


811  TO 


TO  ID 
530  TO  IA 
40  TQ  IB 
523  TO  IC 


C 


ASSIGN 
MM=K 

GO  TO  535 
GO  TO  ID 
ASSIGN  914 
ASSIGN 
ASSIGN 
ASSIGN 
l SL  =- 1 
GO  TO  92  ' 

ISL  =  0 
ALR  = ALRS 
ALI  =  ALIS 
I  T=  1 

__  DO  504  I  =  1  , N 
NEXT  TWO  INSTRUCTIONS  DESIGNATE  INITIAL 
PRECEDING  EIGENVECTOR  SOLUTION 
XR ! I ) =GR ! I i LOOP ) 

XI  {  I)=GI  !  I*£00.P) 

DO  5  1  =  1  » M  a 
AR!  I  ,  I ) =  A R (I , I ) - ALR 
Alt I»I)=AI< I»I J-ALI 
I  J  =  1 


917 

811 


523 

9 


403 


504 

4 


GUESSES  FOR  EIGENVECTOR 


10 


1  1 


12 

13 

106 


1C9 


14 


B  I  G  =  0 . 

DO  13  1=1, N 
YR! I ) =0 . 

YI ! I ) =0 . 

DO  11  J=1,N 

YR! I ) =YR ! I ) + AR { I »J)*XR(J)-AI( I , J ) «  X I  (J) 
YI !  I ) =Y I ( I ) +A I ( I , J ) *XR  t  J )  +  AR { I , J) »XI ( J  ) 
AM= YR ( I)*YR( I )+YI ( I )*YI tl) 

IF  (AM-BIG)  13,13,12 
8 1 G  =  AM 
J  J=  I 

CONTINUE 

IF (BIG)  109,106,109 

I C  T  =  1 000 

JJ  =  N 

ISL  »  1 

GO  TO  92 

RQNR=0 . 

RQN 1=0. 

R0D=0 . 

DO  14  1=1, N 

RQNR=RQNR  +  XR! I ) *YR !  I ) +  X  I { I ) *  Y I ( I) 

RQN  I  =  RQN  I  +  XR  (  I  )  *YI  (  I  5  -XI  (  I  )*YP.  C  I  ) 
RQD=RQD+XR ( I) *XR ( I ) +XI ( I ) - v 
AMUR=RQNR/RQD 


-XI ( I ) 


4-2 


AMU  i  -RQNI  /RGD 
A  t :  M  =  A  M  U  R  -  A  M  U  R  A  M  U I  *  A  M  U  I 
81  TS  =  0„ 

08  15  I  =  1  ,  N 
150TS=TS- 

1  j  y  T  { 

'r.01  l  6'  1  =  1  rN 

XRt  I  )=(  YR(  JJ)  e-YR(  r  )  +  YI  (  J.  •  ..YI  {  n  )/Bi  G 


5  I  =  1  ,  N 

■  S  +  {  V  R.C  I  )  -  AMUR*  XR  (  I  )+AFUl*XI  (  I 
I ) -AMUR* X  I { I ) - AMU  I »X  R ( I } ) **2 

6  I  =  1  f  N 


) )* *2+ 


1  6 

i  1.1 
13 
.19 

20 


fRi  i  n/ciG 


310 

99 

100 

29 

535 


21 


22 


X!  (  I  }  =  { YR{ JJ)*YI ( I ) -V I 
XRC JJ  }  =  1 .0 
XiiJJhC.C 

(TS  /  RQD,-  10.E-4)  20,20,18 
1 1-  (  I J  -  20)  19,20,20 

GO  fn’io 
I C  T= I J 

M I T  2  =  IJ  +  2  0 
AlR=Atf UR+ALR 
All  =  AMU  I  4-AL  1 
'.V-N 

OG  310  1=1 ,N 
ARC  I,  I  )  =  AR  {  I  ,  I  )  -  Af'lJR 
AI(  I,  I)=AI ( I , I )-AMUI 
GO  TO  29 
DO  ICO  I  =  1,N 
ARC  I,  I )  =  AR t I rI ) - ALR 
Aid,  I  )  =AI  (1,1  )  -  A  L I 
i  J=I  J-H 
00  27  1=2, MM 
I M 1  =  1-1 
DO  27  J  =  1  ? I M 1 

FM= A  R  C I , J ) * AR ( I , J  5  +  AI ( I, J ) »AI { I , J ) 

SO" AR { J, J  5* ARC  J, J  5  +  AI t J, J )*AI ( J. J ) 
r  C  F  M- S  M )  24,24,22 
<0  23  K  =  J,MM 
T  1  =  AR  {  J  ,  K  )  •. 

T  2  =  A  I  (  J  ,  K  ) 

ARC  J,K)=AR{ I , K) 

AIC J,K>=AI( I , K ) 

AR (  I , K  )  =  T 1  ' 

23  A  I C I , K) =T2, 

T  1  =  X  R  (  J  ) 

T  2  =  X  I  (  J  ) 

XR £  J ) =  X R (  I  ) 

XI  (  Ji=XI  U) 

XR  { I ) =T 1  ; 

XI { I ) =12  t 
T 1  =  F  M 

f;-’=sm 

s;;=  n 

IF  (SM)  25,27,25 
IF  ( F M )  90  »  27  - °0 

jarc  i,H!  ;  i  J  )  +AIC  I  ,  J  5  *  AI  (  J ,  J  ) )  /SM 
FI =i£Rv J,25* Ai { I , Jl-ARC I, J }*AI ( J, J ) ) /SM 
DO  ^  o  K=J«MM 

ARf I ,KJ=AR( 1 ,K)-RR*AR( J,K)*RI«AI C J,K ) 

(  I  I  1  I  ,  KJ-RRe-A1  {  J  ,  K)  -RI*AR{  J,.<  ) 

AR ( I , j } =0 . 

AT C I , J ) =0* 

XR  C  I ) =XR ( I ) -RR*/R ( J )+RI*XI(J) 

XI {  I 5=XI  ( I )-RR*XI { J)-RI»XR( J) 


24 

25 
90 


26 

27 

550 


75  0 
751 

23 


GO  TO  CA 
SMALL  =  1 000  * 

DO  28  K= 1 , MM 
I  K  K  =  K 

T2=AR(KfKJ.AR(K,K)+AI(KfK>.AI(K,K) 
Ir  (115  roO t 752, 750 
IF  ( T 1 -SMALL  5  75  1,23,28 
SM  ALL  =  T 1 

r:= i< 

CONTINUE 
GO  TO  IB 


ko 


752  12  =  IKK 

IF  USD  753,30,30 

30  157=1 

!l i -2 000 

00  974  1=1, MM 
XR!I ) =0.0 
974  XI (  1  )=0.C 

753  YR(  iZ 5  =  1 .0 
YI (  12  5=0.0 
J  J  =  I  2 

R  1  G=  1 . 0 

IF  ( I Z - 1  *  M )  33,32,33 

32  I Z  2  =  2 

GO  70  95 

33  122=12+: 

DO  31  I  =  I Z  2  ,  MM 
YR( I)=0. 

31  Y II  I ) =0 . 

I ZZ=MM-I Z+2 
1  r  U2-1)  95,49,95 

40  122=1 

41  •:  IG-0. 

95  DO  4  6  I  =  I  2 2 , M M 
11=  MM- I  +  1 
XK  = I  I  + 1 
SR  =  0 . 
si=0u 

IF  (I-I)  42,44,42 

42  DO  43  K=KK,MM 

SR  =  SR+AR ( I  1 ,K)*YR(K)-AI ( I  I , K ) * Y I  (K) 

43  SI  =SH-AR  {  I  I  ,K)*YI  (  K)  +  AI  {  I  I  ,  K)  *YR(K) 

44  T 1 =AR J 1 1 , 1 1) * AR {  1 1  ,  III +A  I  ( I  1 , 1  I ) *A I (  1 1 ,  1 1  1 

YR(  1 1  }  =  {  AR(  l  l  ,  I  I  )  *{  XR{  I  I  3-SR1  +AI  (  I  I,  I  I  )  *(XI  (  I  D-SI  )  )  /T1 
YI { 1 1  }  =  { AR(I I v I  I )*( XI (  I  I  3 -S I >-AI ( I  I,  I  I )* (XR(  II )-SR ) ) /T1 
AM=YR (  I  I  )  »YR { 1 1 ) +YI ( I  I ) « Y I U I) 

IF  (AM-BIG)  46,46,45 

45  JJ=I! 

B I G=AM 

46  CONTINUE  " 

49  DO  47  1=1, MM 

X R ( I )=( YR( JJ>*YR( I ) +YI ( JJ5*YI ( I ) )/8IG 

47  XI U}=(YR(4J)+YI (I)-YI ( JJ)*YR( I ) J/BIG 
xr( j j )  =  i .0; ' 

X I  (  J  J  )  =0 . 0  >,  * 

92  00  601  I = 1 , N 
DO  601  J  =  1  , N 

AR  (  I  ,  J  )  =  B  R  (  I  ,  J  ) 

601  AI(I,J)=BI(I,J) 

116  IF  USD  755,50,60 
755  GO  TO  IC 

50  ALR=0. 

A  L  I  =  0  . 

S  U  M  =  0 . 0 

55  DO  52  1=1, N 
YR( I )=0. 

YI ( I )=0. 

CO  51  K=1,N 

YR (  I  5  =  YR ( I ) +  AR( I ,K) *XR(K  5-AI { I ,K)*XI  ( K ) 

51  YI II) =YI { I ) +AR { I , K ) *  X I (  K  )  +  A  I (  I ,  K  )  *  XR  (  K  ) 

ALR=ALR+XR( I  )*YR( I )  +  X  I (  I  )  *Y I ( I ) 

AL I  =  ALI  +XR(  1) «YI ( I ) -X  I ( I  )*YR( I ) 

52  SUM=SUM+XR( I)«XR( I)+XI (I )«XI (I) 

ALR=ALR/ SUM 

AU  =AL  1/  SUM 
AM=ALR«ALR+ALI*ALI 
B3  TS=0 , 

DO  53  1=1 , M 

T I  =  V R  (  I)-ALR*XR(  l)+ALI*XI{  I  ) 

T2  =  Y  I  (  I  )-Al.R*XI  (  I  )-ALI*XR(  I  ) 

53  TS=TS+T1*T1+T2*T2 

93  IF  (IS  /  SUM  -  10.E-14)  60,60,301 
301  IF  (IJ-MI72)  99,400,400 

400  IF  (IT  -  3)  402,990,402 
402  ALR  =  ALR  +  .1 


60 

63 


All  =  All 
I  T  =  I  T  +  1 
CO  10  4 
I  S  L  =  0 


61 


911 

914 

915 

703 


VALUES { 1 , LOOP)  =  ALR 
VALUES { 2 ,  LOOP  )  =  ALI 
IF  UJ-N)  61,65,61 
ti=xr(jj) 

T  2  =  X  I { JJ  ) 

XR [ JJ  )  =  X  R  (  N  ) 

X I  f  J  J  )  =x  I  {  N!) 

X  R  (  N  )  =  T 1 
XI (N)=T2 
CO  63  K=1,N 
1  1  =AR  { 

T  2  =  A  I  {  JJ,  |<  i 
AR{ JJ ,K ) =AR( M,K ) 

A  I  (  J  J  ,  K  )  =A  I  {  N ,  K ) 

AR  ;  N ,  K )  =  T 1 
63  A  I  f.\',KJ«72 
00  62  K  =  1  , N 
1  3  —  A  R  {  K  ,  J  J  ) 

!  2  =  A  I  (  K  ,  ‘J  J  ) 

AR ( K , J  J ) =  AR ( K , N ) 

A I  K, JJ)=AI (K,N) 

AR  {  K ,  N )  =  T1 
A  I  I  K ,  N  }  =  T2 
N=N-  I 

DO  66  I  =  1  ,  M 
DO  66  J  =  1  ,  N 

=  I  ?  J i)-XR{  I  )*AR  (N+l  ,J)+XI(I  }*AIIN+1  n 

DO  600  i  =  J)~XR(I,*AI  (N+l,J)-xni  )*AR(N+llj) 
DO  600  J  =  1  ,  Xi 
BR {  I , J )  =  AR{ I , J ) 
dlf I, J)=AI ( I , J) 

DO  702  1  =  1  M 
DO  T02  J  =  1 1 M 
ARC  I,  J)=CRU,  J) 

AI!  ,Jl=Cin,  J) 

Ir  Il-J)  762,701,702 
AR}  I»  U=ARfcJ,,  r  ) -ALR 

CONliNuI4  '  -  *  1 }“ALI 

ASSIGN  911  'TO 
GO  TO  535 
ASSIGN  530  TO 
I SL=- 1 

00  703  I = I , M 
XR ( I )=0. 

XI { I >=0. 

ASSIGN  753  TO 
ASSIGN  7G4  TO 
GO  TO  530 
ASSIGN  525  TO  IC 
GO  TO  92 


62 

65 


66 


600 

700 


701 

702 


IA 
I  A 


IB 

IC 


916 

704  DO  920  * !  =  i,m 

EIGENVECTORS ’xR^Jnd^X!  FR0M  rRAK'SF0RMATI  ON  MATRICES  GR  AND  GI  FROM 
GRUjLOOP)  =  XR  (  I  } 

Gilt, LOOP }  =  X[(l) 

ASSIGN  40  TO  IB 
LOOP  =  LOOP  +  1 
IfCN-l)  921,67,523 
ALR=AR (1,1) 

AL I =A  I  ( 1  ,  1  ) 

VALUES ( 1 , LOOP )  =  ALR 
nVALUES(2,LC0P)  =  ALI 

GO  TO  700 
921  GO  TO  4000 


920 


525 

67 


433 


XMAX  =  -10.  «*  2  0 
03  1 000  I  =  J,M 

rn  (XMAX  -  VALUE  S ( 1 , I ) )  3,3,1000 
3  X  AX  -  VALUES!  I,  15 

1  <  X  —  I 

1CC0  COM! I  VUE 

IF  (MX-J)  1,4000,1 
i  VALUES (1, MX)  =  VALUES! 1 , J ) 

V A L  U E  S ! 1 , J  5  =  XMAX 
SAVE  =  VALUES (2, MX) 

VALUE S (2, MX) =  VALUES ( 2, J) 

VALUES ( 2, J 5  =  SAVE 
03  3000  K=1,M 
SAVE  =  OR (K, MX  5 
GfUK,MX3  =  GR(KtJ) 

CR { K , J )  =  SAVE 
SAVE  =  G  I  ( K , M X ) 

6I(K,HX)=  6I(K,J) 

3000  GI(K,J)  =  SAVE 
4000  CONTINUE 
GO  TO  992 

990  PRINT  991 

991  FORMAT ( I 5H  NO  CONVERGENCE) 

992  END 

SUBROUTINE  GAUSS3 ( N , E?,A, X,KER> 
DIMENSION  A ( 40 , 40  5 ,  X!40,40) 

DO  1  1  =  1,  N* 

DO  1  J  = 1 , N 

1  X(!tJ)=0„Q 
DO  2  is  =  1  ,  N 

2  x  { :< ,  ;< )  =  U  0 

10  00  34  L  = 1 9 N 
KP  =  G 

2  =  0.0 

DO  1 2  K=  L , N 

IF(Z-ABSF(A(K,L)  )) 1 1, 12,  12 

11  Z  =A8SF  !  A  !  Kj-L  )  5 

j/  p  —  !^  '  r  ■ 

12  CONTINUE  > : 

IF(L-KP) 1  3 »  2  0  }  2  0 

13  CO  14  J=L,N 

Z  =  A ( L , J  5  " 

A ( L , J ) =A ( KP , J ) 

14  A  !  K  P  ,  J  )  =  Z 
00  15  J= 1 »  N 
Z  =  X  (  L  j  J  ) 

X ! L , J  3  =X i KP »  J ) 

15  X ! KP  f  J } =Z 

20  IFtASSF! A!L,L ) ) -EP ) 50 , 50 , 30 

30  IF(L-N) 51,34,34 

31  L  P  1  =  L  -J- 1 

DO  36  K=  LP 1 , N 
I F  {  A  {  K ,  L  ) ) 32,36,32 

32  RATIQ=A(K,L)/A(L,L) 

DO  33  J  =  l. P 1  ,  N 

33  A(K, J)=A(K, j)-RATIO*A(L» J) 

DO  35  J  = I , N 

3  5  X(K, J)=XvK»J)-RATIO»X(L,  J) 

36  CONTINUE 

34  CONTINUE 

40  Du  43  1  =  1  ,N 
I  I  =  N  +1-1 

DO  43  J  =  I , N 
S  =  0.0 

IF!  1 1  - M } 4 1 , 43,43 

41  I  I P 1  =  1 1  +  1 

DO  42  K=  1 1 P 1 , N 

42  S  =  S  +  A  { II , K ) *  X ( K , J ) 

43  X(H»J5  =  (X(II,J)-S)/A<II,II) 

KER=  1 

GO  TO  15 
50  KER=2 

‘75  CONTINUE 
END 
END 

430 
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GRAPHS 
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48 


•  49 


£0 


V2 


ap 
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